{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--BOOK_INFORMATION-->\n",
    "<img align=\"left\" style=\"padding-right:10px;\" src=\"figures/PHydro-cover-small.png\">\n",
    "*This is the Jupyter notebook version of the [Python in Hydrology](http://www.greenteapress.com/pythonhydro/pythonhydro.html) by Sat Kumar Tomer.*\n",
    "*Source code is available at [code.google.com](https://code.google.com/archive/p/python-in-hydrology/source).*\n",
    "\n",
    "*The book is available under the [GNU Free Documentation License](http://www.gnu.org/copyleft/fdl.html). If you have comments, corrections or suggestions, please send email to satkumartomer@gmail.com.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [NDVI](06.08-NDVI.ipynb) | [Contents](Index.ipynb) | [Unsupervised Classiﬁcation](06.10-Unsupervised-Classiﬁcation.ipynb) >"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.9 调用GDAL作为外部命令\n",
    "\n",
    "GDAL库提供了许多有用的程序。例如，修改GIS数据的坐标系统、分辨率和范围等等。完成的列表可以在`http://www.gdal.org/gdal_utilities.html`找到。在本节中，我们将改变数据的分辨率、范围和坐标系统。不幸的是，这些程序没有可用的Python绑定，所以我们不能像使用其他Python程序那样使用它们。这些GDAL程序在命令行模式下工作，Python提供了一种从Python中运行外部命令的方式。所以，我们将使用这个来利用GDAL库的强大功能。在跳入GDAL之前，我们将简要地看一下如何从Python中运行外部的命令。在Linux上，我们通过使用命令`ls`来得到一个目录中文件列表，`ls -al`也提供了这些文件/目录的细节。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "from subprocess import call\n",
    "dir_list = call([\"ls\", \"-l\"])\n",
    "print(dir_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "请注意，我们将该参数作为列表的一个单独成员传递给该命令。随着参数数量的增加，这会变得很混乱。在linux shell上，我们简单地写了`ls -al`，为什么我们不能采用相同的方法在这里做？让我们试着像这样写。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'ls -al'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-faadbd7dc60f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdir_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"ls -al\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/envs/datascience/lib/python3.5/subprocess.py\u001b[0m in \u001b[0;36mcall\u001b[0;34m(timeout, *popenargs, **kwargs)\u001b[0m\n\u001b[1;32m    245\u001b[0m     \u001b[0mretcode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"ls\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"-l\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    246\u001b[0m     \"\"\"\n\u001b[0;32m--> 247\u001b[0;31m     \u001b[0;32mwith\u001b[0m \u001b[0mPopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpopenargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    248\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    249\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/datascience/lib/python3.5/subprocess.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds)\u001b[0m\n\u001b[1;32m    674\u001b[0m                                 \u001b[0mc2pread\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc2pwrite\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    675\u001b[0m                                 \u001b[0merrread\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merrwrite\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 676\u001b[0;31m                                 restore_signals, start_new_session)\n\u001b[0m\u001b[1;32m    677\u001b[0m         \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    678\u001b[0m             \u001b[0;31m# Cleanup if the child failed starting.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/datascience/lib/python3.5/subprocess.py\u001b[0m in \u001b[0;36m_execute_child\u001b[0;34m(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)\u001b[0m\n\u001b[1;32m   1287\u001b[0m                             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1288\u001b[0m                                 \u001b[0merr_msg\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m': '\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mrepr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0morig_executable\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1289\u001b[0;31m                     \u001b[0;32mraise\u001b[0m \u001b[0mchild_exception_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merrno_num\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merr_msg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1290\u001b[0m                 \u001b[0;32mraise\u001b[0m \u001b[0mchild_exception_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr_msg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1291\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'ls -al'"
     ]
    }
   ],
   "source": [
    "dir_list = call(\"ls -al\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python以这种方式报错。要做到这一点，我们应该告诉Python，像在shell中一样，通过提供参数`shell=True`来完成输入。看看下面的例子。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "dir_list = call(\"ls -al\",shell=True)\n",
    "print(dir_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里给出了想要的输出。当你需要将GDAL作为外部命令行调用时，这足以了解`subprocess`模块使用外部命令行。让我们先改变数据的UTM投影，然后从影像中选择子窗口，并最终改变数据的分辨率。让我们先看看如何在命令模式下执行此操作。在Python脚本中使用此命令之前，首先在shell中尝试您的命令始终时一种好的做法。我们将使用ASTER DEM来操作。ASTER数据在地理坐标中提供。我将使用`ASTGTM_N11E076.tif`影像，这个印度南部kabini流域的DEM。该文件的路径在我的文件中是`/home/tomer/berambadi/DEM/ASTGTM_N11E076_dem.tif`,将其替换成您的文件路径。首先我们将看到这个tif文件的路径。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "$ gdalinfo /home/tomer/berambadi/DEM/ASTGTM_N11E076_dem.tif\n",
    "Warning 1: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order\n",
    "Driver: GTiff/GeoTIFF\n",
    "Files: /home/tomer/berambadi/DEM/ASTGTM_N11E076_dem.tif\n",
    "Size is 3601, 3601\n",
    "Coordinate System is:\n",
    "GEOGCS[\"WGS 84\",\n",
    "    DATUM[\"WGS_1984\",\n",
    "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
    "            AUTHORITY[\"EPSG\",\"7030\"]],\n",
    "        AUTHORITY[\"EPSG\",\"6326\"]],\n",
    "    PRIMEM[\"Greenwich\",0],\n",
    "    UNIT[\"degree\",0.0174532925199433],\n",
    "    AUTHORITY[\"EPSG\",\"4326\"]]\n",
    "Origin = (75.999861111111116,12.000138888888889)\n",
    "Pixel Size = (0.000277777777778,-0.000277777777778)\n",
    "Metadata:\n",
    "    TIFFTAG_DOCUMENTNAME=created at\n",
    "    TIFFTAG_IMAGEDESCRIPTION=SILC TIFF\n",
    "    TIFFTAG_SOFTWARE=IDL 6.3, ResearchSystems,Inc.\n",
    "    TIFFTAG_DATETIME=2008:10:28 23:54:25\n",
    "    TIFFTAG_XRESOLUTION=100\n",
    "    TIFFTAG_YRESOLUTION=100\n",
    "    TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)\n",
    "    AREA_OR_POINT=Area\n",
    "Image Structure Metadata:\n",
    "    INTERLEAVE=BAND\n",
    "Corner Coordinates:\n",
    "Upper Left ( 75.9998611, 12.0001389)(75d59'59.50\"E,12d 0' 0.50\"N)\n",
    "Lower Left ( 75.9998611, 10.9998611) ( 75d59'59.50\"E,10d59'59.50\"N)\n",
    "Upper Right ( 77.0001389, 12.0001389) ( 77d 0' 0.50\"E,12d 0' 0.50\"N)\n",
    "Lower Right ( 77.0001389, 10.9998611) ( 77d 0' 0.50\"E,10d59'59.50\"N)\n",
    "Center( 76.5000000, 11.5000000) ( 76d30' 0.00\"E,11d30' 0.00\"N)\n",
    "Band 1 Block=3601x1 Type=Int16,ColorInterp=Gray\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在`gdalinfo`输出的大多数参数都是自我解释的。我们将分析所需要的东西，例如坐标系、分辨率和角坐标。看到线`AUTHORITY`，值为`EPSG,4326`。网站`http://spatialreference.org/`可以看到输出中与地理参考有关的所有术语的含义。对于`EPSG:4326`,网页是`http://spatialreference.org/ref/epsg/4326/`，它提供了关于这个投影的所有细节。如果影像的投影系统被很好的定义，那么在投影到一个新系统中的投影就无关紧要的。数据的分辨率是0.000277度，数据的范围在最后一行的输出中给出。\n",
    "\n",
    "改变投影系统的命令是："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "$ gdalwarp -r lanczos -t_srs '+proj=utm +zone=43 +datum=WGS84' \\\n",
    "> /home/tomer/berambadi/DEM/ASTGTM_N11E076_dem.tif /home/tomer/berambadi/DEM/temp0.tif\n",
    "Warning 1: TIFFReadDirectoryCheckOrder:Invalid TIFF directory;\n",
    "tags are not sorted in ascending order\n",
    "Processing input file /home/tomer/berambadi/DEM/ASTGTM_N11E076_dem.tif.\n",
    "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "您可以通过再次发布命令`gdalinfo`，来查看新创建的文件`temp0.tif`的投影系统信息。现在让我们从这些数据中，改变坐标（xmin = 664000，YMAX = 1309000，XMAX = 685000，Ymin = 1294000）的子窗口，命令是："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "$ gdal_translate -a_ullr 664000 1309000 685000 1294000 \\\n",
    "> -projwin 664000 1309000 685000 1294000 \\\n",
    "> /home/tomer/berambadi/DEM/temp0.tif /home/tomer/berambadi/DEM/temp1.tif\n",
    "Input file size is 3595, 3645\n",
    "Computed -srcwin 1807 601 688 492 from projected window.\n",
    "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在你又可以使用`gdalinfo`来查看新创建的`tem1.tif`的文件的信息。你会看到数据的分辨率大约是30m。现在让它达到50m。从这些数据中，改变坐标(xmin=664000, ymax=1309000, xmax=685000, ymin=1294000)的子窗口。这样做的命令是`gdalwarp`，这里有两个必要的参数`:r`和`-tr`。`r`参数定义了插值算法，`tr`定义了水平和垂直方向上的分辨率。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "$ gdalwarp -r lanczos -tr 50 50 \\\n",
    "> /home/tomer/berambadi/DEM/temp1.tif \\\n",
    "> /home/tomer/berambadi/DEM/dem_50m.tif\n",
    "Creating output file that is 210P x 150L.\n",
    "Processing input file /home/tomer/berambadi/DEM/temp1.tif.\n",
    "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在通过使用`gdalinfo`命令给出关于这个最终的数据集的信息如下："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "$ gdalinfo /home/tomer/berambadi/DEM/dem_50m.tif\n",
    "Driver: GTiff/GeoTIFF\n",
    "Files: /home/tomer/berambadi/DEM/dem_50m.tif\n",
    "Size is 420, 300\n",
    "Coordinate System is:\n",
    "PROJCS[\"WGS 84 / UTM zone 43N\",\n",
    "GEOGCS[\"WGS 84\",\n",
    "DATUM[\"WGS_1984\",\n",
    "SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
    "AUTHORITY[\"EPSG\",\"7030\"]],\n",
    "AUTHORITY[\"EPSG\",\"6326\"]],\n",
    "PRIMEM[\"Greenwich\",0],\n",
    "UNIT[\"degree\",0.0174532925199433],\n",
    "AUTHORITY[\"EPSG\",\"4326\"]],\n",
    "PROJECTION[\"Transverse_Mercator\"],\n",
    "PARAMETER[\"latitude_of_origin\",0],\n",
    "PARAMETER[\"central_meridian\",75],\n",
    "PARAMETER[\"scale_factor\",0.9996],\n",
    "PARAMETER[\"false_easting\",500000],\n",
    "PARAMETER[\"false_northing\",0],\n",
    "UNIT[\"metre\",1,\n",
    "AUTHORITY[\"EPSG\",\"9001\"]],\n",
    "AUTHORITY[\"EPSG\",\"32643\"]]\n",
    "Origin = (664000.000000000000000,1309000.000000000000000)\n",
    "Pixel Size = (50.000000000000000,-50.000000000000000)\n",
    "Metadata:\n",
    "AREA_OR_POINT=Area\n",
    "Image Structure Metadata:\n",
    "INTERLEAVE=BAND\n",
    "Corner Coordinates:\n",
    "Upper Left ( 664000.000, 1309000.000) ( 76d30'19.72\"E, 11d50'14.13\"N)\n",
    "Lower Left ( 664000.000, 1294000.000) ( 76d30'17.06\"E, 11d42' 5.94\"N)\n",
    "Upper Right ( 685000.000, 1309000.000) ( 76d41'53.51\"E, 11d50'10.21\"N)\n",
    "Lower Right ( 685000.000, 1294000.000) ( 76d41'50.51\"E, 11d42' 2.07\"N)\n",
    "Center\n",
    "( 674500.000, 1301500.000) ( 76d36' 5.21\"E, 11d46' 8.15\"N)\n",
    "Band 1 Block=420x9 Type=Int16, ColorInterp=Gray\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个最终的数据具有所需的UTM坐标、范围和分辨率。\n",
    "\n",
    "让我们在Python中完成所有这些步骤。首先我们将制作一个字符串，它基本上是shell中使用的命令。然后我们将使用`call`函数来运行它。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# specify file names\n",
    "input_file = '/home/tomer/berambadi/DEM/ASTGTM_N11E076_dem.tif'\n",
    "temp_file0 = '/home/tomer/berambadi/DEM/temp0.tif'\n",
    "temp_file1 = '/home/tomer/berambadi/DEM/temp1.tif'\n",
    "out_file = '/home/tomer/berambadi/DEM/dem_50.tif'\n",
    "# convert from geographical co-ordinates to UTM zone 43\n",
    "shell_command= \"gdalwarp -r lanczos -t_srs '+proj=utm +zone=43 +datum=WGS84' %s %s\"%(\n",
    "input_file,temp_file0)\n",
    "returncode = call(shell_command, shell=True)\n",
    "berambadi = ' 664000 1309000 685000 1294000 '\n",
    "# cut the area around Berambadi watershed\n",
    "shell_command = \"gdal_translate -a_ullr %s -projwin %s %s %s\"%(berambadi,\n",
    "berambadi, temp_file0, temp_file1)\n",
    "returncode = call(SC, shell=True)\n",
    "# changing the resolution\n",
    "res = 100\n",
    "shell_command= \"gdalwarp -r lanczos -tr %s %s %s %s\"%(res, res,\n",
    "temp_file1, out_file)\n",
    "returncode = call(SC, shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GDAD中的某些程序并不覆盖现有的文件，或者如果文件存在相同的名称，可能会创建一些有趣的数据集。因此，最好确保没有用于创建新数据的现有文件。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
